Chapter 5 Areal Unit Data
This section will focus on mapping areal unit variable, where we have a variable of interest that we want to display for each region of interest. For this section we will look at building a map of SIMD quintile across the Lothian Health Board. In order to do this we will need to:
Read in the spatial information using the function
read_sf()from the{sf}packageRead in the data that we want to plot using the
{readr}package (or a relevant package for reading in your data file)Join with a standard data set, and filter using functions from the
{dplyr}packageProduce maps using either the
{ggplot2}package, or the{leaflet}package
We’ll start by loading all of these packages - remember if you are only using {ggplot2} you don’t need to load {leaflet} and similarly if you are only using {leaflet} you don’t need to load {ggplot2}.
library(sf) # reading in the spatial information
library(readr) # reading in data (this may change depending on file type)
library(dplyr) # manipulating the data
library(ggplot2) # plotting a static map
library(leaflet) # plotting an interactive map5.1 Reading in Shapefiles
The spatial information about the regions that you are trying to plot is contained within the shape files. This is a collection of files all with the same name and different file extensions that contain information on the regions, and how they should be plotted. For R you only need three of these files with the extensions .shp, .dbf, and .shx. Each of these serves a particular function:
The
.shpfile (shape file) contains information on the boundaries of the regions. Specifically this is a collection of sets of coordinates that make up the boundary.The
.dbffile (database file) contains data that relate to the boundaries. This can include multiple different variables, but at a minimum must contain a column with region IDs. Usually this file will also contain information on things like the region areas, though these typically are not used.The
.shxfile (shape index file) contains indexing information that allows the information in the shape file to be linked to the data in the database file. It also has implications on plotting order (though this usually isn’t particularly important).
Since these files are all linked, they are read in together using the function read_sf(), where the only argument is the file path for the .shp file. The code below reads in information on data zone boundaries and names the resulting object spat:
path <- "/conf/linkage/output/lookups/Unicode/Geography/Shapefiles/Data Zones 2011/"
spat <- read_sf(paste0(path, "SG_DataZone_Bdry_2011.shp"))When you read this in, you should notice that it appears to be exactly the same as any other data frame in R. It pretty much is, which means that you can join, and filter and manipulate the data in spat in exactly the same way as any other data frame. The thing that makes spat spatial is that is has a special column called GEOMETRY which contains objects that are POLYGON or MULTIPOLYGON types defining the region associated with each observation.
5.2 Joining with data
The next step is to read in the SIMD information - this step is exactly the same as it normally would be. The information in this case is stored in a .rds file, so we will make use of the function read_rds() from the readr package. The column Datazone2011 is being renamed to DataZone to match the column name in the spatial information.
data_path <- "/conf/linkage/output/lookups/Unicode/Deprivation/"
SIMD <- read_rds(paste0(data_path, "DataZone2011_domain_level_simd.rds")) %>%
rename(DataZone = Datazone2011)Now that we have all of the information required to map, we need to join the two data sets. We do this in the usual way using one of the join functions from {dplyr}. One important note here is that the object containing the spatial information must be the first argument in your join function, if it isn’t your resulting data will not have spatial characteristics:
class(left_join(SIMD, spat))## [1] "tbl_df" "tbl" "data.frame"class(left_join(spat, SIMD))## [1] "sf" "tbl_df" "tbl" "data.frame"You can see that the class of the combination of SIMD and spat does not have the characteristic sf while the class of the combination spat and SIMD retains this characteristic.
Since we are only interested in plotting over regions in the Lothian health board, we should also filter these out. We can do this as usual using the filter() function from {dplyr}. So obtaining our joined and filtered data set as follows:
SIMD_map <- spat %>%
left_join(SIMD) %>%
filter(hb2019name == "NHS Lothian")Now that we have the data joined and filtered appropriately, we can start producing maps.
5.3 {ggplot2}
To produce a map in ggplot2 with any spatial data (including point data) we use the geometry function geom_sf(), within which we can specify aesthetics like fill and colour to represent observations. So, to map SIMD quintile, which is represented by the column simd202_11:
ggplot(SIMD_map) +
geom_sf(aes(fill = simd202_11))
This is fine, but we can customise it a bit more by changing the colour scale, and naming the fill scale more appropriately. To help to visualise the small regions, I’ll also set the border colour to be set according to SIMD. Maps also usually look nicer when they are presented using theme_void() (although this is up to your own personal preference!):
ggplot(SIMD_map) +
geom_sf(aes(fill = simd202_11, colour = simd202_11)) +
scale_fill_viridis_c("SIMD Quintile", option = "C", aesthetics = c("fill", "colour")) +
theme_void()
Although this map is static, we can use ggplotly() to introduce a small amount of interactivity. The column has been renamed to improve the labels, and tooltip has been specified to reduce duplication. There is also a dummy aesthetic DZ which allows us to show the data zone name:
p <- SIMD_map %>%
rename(`SIMD Quintile` = simd202_11) %>%
ggplot() +
geom_sf(aes(fill = `SIMD Quintile`, colour = `SIMD Quintile`, DZ = DataZone)) +
scale_fill_viridis_c("SIMD Quintile", option = "C", aesthetics = c("fill", "colour")) +
theme_void()
plotly::ggplotly(p, tooltip = c("DZ", "fill"))